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CHAPTER 1 

Introduction - Qualitative Review 
1.1 Motivation 

Interest in the strength of naturally occuring ice and 
the ultimate size to which marine structures must he designed 
has become a major topic as the last boundaries of easily 
accessible petroleum resources are used. The extension of 
engineering practice into the world's colder regions is 
subsequently an effort that is intimately related to this 
perceived need. 

To satisfy a curiosity regarding this area of offshore 
design, the author set forth to discover a rational approach 
to establishing a meaningful standard. It was hoped that 
an appropriate theory, perhaps based on statistical patterns, 
was available that was similar to the many comprehensive 
theories regarding wind, wave, and seismic loading. Unfor- 
tunately, no comprehensive or universally accepted theory 
was found. However, many prominant investigators, including 
Assur [1 ][2][3]* Bercha [4][5][6]* Hirayama et al. [14] 
[15]* and Korzhavin [19] [20]* have proposed numerous models 
and approaches to the calculation of ice strength and the 
quantitative prediction of necessary design size. Some of 
these models will be reviewed in subsequent chapters. 

Consequently, the motivation for this thesis is to 
recognize and segregate an accurate and useful theory for 
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the prediction of ice pressures. 

1.2 Scope 

The purpose of this investigation is to analyze the 
ice-structural interaction models presently proposed for the 
case of a vertical circular pile. The circular pile was 
chosen primarily because it is the most useful geometric 
form in areas of high tidal variance. Other authors , 
notably Tryde [ 30 ] and Bercha [ 5 ]* have examined the other 
promising forms , the inclined wedge and the cone. 

The investigation will take the form of the following 
steps : 

1) Qualitatively compare the available analytical 
models . 

2) Either improve upon the most promising model or 
construct a new model. 

3) Analyze sensitivity to some important parameters. 

4 ) Analyze some important ice properties affecting 
the circular pile interaction problem. 

1.3 Structural Failure 

The earliest designs contemplating an ice environment 
were highway and railroad bridges in North America and the 
Soviet Union. Some authors, notably Korzhavin [ 19 ] and 
Watts [ 32 ], have summarized structural collapses known to 
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"be due to ice forces. These failures have a common thread; 
they were all caused by unusual ice jams which, coupled 
with an early thaw, produced unplanned loadings. Insights 
gained from a careful reading of these descriptions indicate : 

- very few failures actually occurred. This evidence 
implies that early engineers had either very accurate 
analytical skills (hardly likely) or that early 
designs were overly conservative. 

- some failures were by gradual ice-erosion. 

- when failure did occur, it was due to an unusual 
meteorological complication or no allowance for 
ice at all. 

No major marine structural failures, apart from 
ice-bound ships, were found in the literature. Perhaps this is 
because this hostile environment has yet to be seriously 
challenged. Ice damage on coastline facilities, breakwaters, 
docking facilities, and lights, have been recorded. [26] 

1.4 Application 

Applications of ice engineering extend not only to the 
offshore industry but 'over a wide range of related problems. 
The knowledge gained through the solution of marine ice 
uncertainty can apply to the following range of actual 
proposals ; 

- over ice transportation 

- ice airfields 

- submarine surfacing through ice 
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- icebreaker technology - navigation 

- ice island construction for use as a moblie platform 

- iceberg towing 

The offshore oil industry has perhaps the greatest 
motivation in subduing arctic conditions. Although 
experience to date is not extensive, the few structures 
presently deployed emphasize the continuing need for more 
accurate analysis. 

Present North American interest is concentrated in 
Cook Inlet, Alaska, [26], the North Slope, and the Beaufort 
Sea [9]. Arctic conditions are most severe, as it will 
be shown in the next few sections, and there are designs 
under construction which do not even consider a passive 
survival of a structure. These proposals include a vibra- 
tory icebreaker motion [7], ice melting by internal heating 
[17], air cushion mobile platform [22], construction of 
protective soil berms [9], or even mobile platforms that 
only drill in optimum conditions. 

Scandinavian, Soviet, and Argentinian oil interests 
are not represented in the literature. 

1.5 Historical Treatment 

Analytical treatment of ice pressures and winter loading 
was first presented in fresh water civil works, especially 
dams and bridge piers. Early investigators faced the problems 
by correctly deducing that no matter what the loading condi- 
tion or geometric possibilities, ice strength could not exceed 
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its crushing strength. Thus a limiting value of O cr was 
imposed into the following relationship: 

f*MAX ~bh 6>CJ* (1«1) 

b represents pier width and 6 CR crush strength of ice 
is force exerted on structure 

A special study of the A.S.C.E. in 1931 (Committee on 
Power Division) recommended a crush strength of 400 psi [32] • 
A review of Russian, Polish, Canadian, and U.S. military 
design codes was conducted by Watts [32], and it reveals 
that most design criteria are simply a refined empirical 
version of (1.1). 

Design criteria for offshore structures in the much 
more quantitatively unknown sea environment do not simplis- 
tically prescribe this relation. The paragraph relating to 
design for sea ice forces for fixed offshore structures from 
Pet Norske Veritas 197^ Rules states [10] : 

B^Ol. For structures intended to be installed in areas 
where ice hazards may exist, relevant statistical data 
for the area in question are to be submitted. The ice 
conditions are to be described with particular attention 
to: 

- concentration and distribution of ice 

- types of ice (ice floes, ice ridges, rafted ice, etc.) 

- mechanical properties of ice 



14 



- mean thickness of ice floes 

- average drifting speed and direction of ice 

- probability of encountering icebergs 

- tidal range 

This rather nebulous requirement is in stark contrast 
to the curves, tables, formulae, and data presented for other 
force loadings. It is here that the lack of analysis in 
the state of the art is most keenly felt. 

To increase the body of knowledge presently available 
in the science regarding ice pressure loadings, numerous 
investigators in the last decade have proposed many alternate 
procedures. The more promising approaches will be analyzed 
in a subsequent chapter. 

The proposed solutions to the ice pressure problem 
fall basically into four broad categories: 

1) Classical elastic thin plate theory 

2) Finite element or finite difference numerical solutions 

3) Empirical formulation 

4) Model basin testing 

Each approach has merits as well as severe, restrictive 
limitations. For instance, the classical theory assumes 
thin plate elastic behavior (plane stress or plane strain, 
isotropic, homogeneous, elastic) which is most certainly not 
the complete behavior of sea ice. Yet the solution can yield 
a fair degree of accuracy within certain limiting conditions 
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and provide a closed form expression which is easy to use 
and easily reproducible. 

1.6 Importance of Ice Forces 

To insure that the reader is qualitatively aware of the 
enormous magnitude that ice pressures can exert, a simple 
model is presented in Figure 1-1. 

A monopod platform with a radius at ice of 8 feet is 
subjected to an ice load driven by wind or current shear. If 
the extent of the floe is great enough, it is obvious that 
forces well in excess of local ice failure can be generated. 
Using a of 200 psi and formula (1.1), adjusted for 

projected frontal area, the ice force can be: 

Fmax ~ Gcftbh 0*90 "tons 

This approximately 700 tons is concentrated on a three 
foot slice of the leg, presenting severe structural consid- 
erations. If this structure were located in 80 feet of 
water, the overturning moment could be 20,700 foot-tons, a 
serious concern to foundation and soil resistance as well as 
flexural strength of the leg. 

Perhaps an even more serious consideration which 
has recently been addressed [21 ] [25] [26] Is the fact 
that this force can thrust at a frequency close to one hertz. 
Apart from resonant considerations, the fatigue strengths 
of the cold metal and cold weldments become of serious 



concern 



16 




FIGURE 1-1 



Ice interaction with monopod type structure 
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Peyton [26] has summarized some important requirements 
for the ice force design of marine structures. These include 

- Horizontal forces due to crushing, bending, buckling, 
or shearing of the ice sheet (quasi-static) 

- Impact forces (dynamic) 

- Weight of attached sea ice at low tide 

- Buoyant lift of above at high tide 

- Vertical components of sheet failure 

- Diaphragm bending during water level change 

- Ice accretion by spray 

- Thermal expansion of ice in joints 

- Rubble in structural framing carried by ice 

- Abrasion by moving ice 

It can be seen that in a cold environment, ice forces 
can be large and obvious as well as an insidious problem. 

In the next section the enormity of the design equation 
and the limits of present knowledge are presented. 

1.7 Ice and the Enormity of the Design Problem 

Ice is a naturally occurring material which forms in a 
convenient plate shape located at the surface of the liquid- 
air interface. Ice formed on a small fresh water pond can 
be clear, homogeneous, constant in thickness and have a sen- 
sible temperature distribution. In this case the engineer 
could possibly apply some of the more well-known classical 
plate solutions to find the ice behavior and predict with a 
fair degree of confidence any desired solution to his prob- 
lem. 
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Sea ice is not that well behaved, however. Consider 
the situation shown in Figure 1-2 which depicts a profile 
of arctic coastal ice. A bottom founded structure could 
be located anywhere in the picture. 

Figure 1-3 presents a careful summary of the unknown 
parameters affecting the design equation. A high degree 
of interdependence is characteristic. 

Figure 1-4 presents four logical paths to the final 
design. All have been used in successful designs. Again, 
each path has its advantages and disadvantages. 

1) From known environmental conditions - this is 
perhaps the ideal situation; a successful design 
could possibly be created anywhere at minimal 
cost, once perfected. Of course, the lack of 
knowledge regarding the environment at each locale, 
the rudimentary state of analytical tools, and the 
lack of correlation with known quantities makes 
this procedure years away, and perhaps unattainable. 

2) In situ - This will provide better data, but still 
suffers from analytical model imperfection, as above. 
In addition, the gathering of this data can be 
costly. 

3) Near full scale - a very practical approach, but 
it too can be costly, nearly as much as the final 
product in some cases [25] • It does offer the 
advantage of risking a small, inconsequential failure 
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FIGURE 1-2 

Typical winter features, Beaufort Sea 
(after Croasdale [9j ) 
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FIGURE 1-4 

Logical procedure for 

design criterion determination 
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over a larger* more serious accident. 

4 ) Previous experience - This approach would yield 
the highest confidence level* but, unfortunately, 
little has been done and some of what has is 
proprietary information. Here again an overly 
conservative design would propagate other overly 
conservative designs. 

Without presenting a complete review of the parameters 
in Figure 1-3 or the literature on Figure 1 - 4 * both of which 
are beyond the scope of this thesis* it is hoped that an 
appreciation of the enormity of the design effort involving 
ice can be gained. 
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CHAPTER 2 

Summary of Some Sea Ice Material Properties 



2.1 Purpose of this Chapter 

The purpose of this chapter is to present a brief 
outline of some important sea ice parameters and properties 
for use in developing the analytical model presented later 
in this thesis. For a thorough treatment of ice properties , 
the interested reader is referred to references [12] , [13] , 
[18], and [33]. 



2.2 Problem Formulation 

Sea ice can be considered to be sensitive to temperature, 
salinity, crystal orientation, and seasonal variation. It 
can be considered to behave elastically for a short load 
duration and a moderate load magnitude. Katona and Vaudrey 
[18] describe an approach to an appropriate field theory 
which is suff iciently general to apply to the case of a 
circular pile - ice sheet interaction problem. The basic 
parameters are presented in Figure 2-1. 

Katona and Vaudrey further reduce the list of material 
parameters (category I of Figure 2-1) to a functional form 
which contains only two parameters : 

t 



where 



<V 



O 



( 2 . 1 ) 



Oj: Is the stress tensor at spatial point X 

J 

and time "t 
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FIGURE 2-1 

Parameters to define an appropriate field theory 
(after Katona and Vaudrey [j 8j ) 
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rjj k | is a symmetric fourth order tensor functional 
allowing complete anisotropic formulation 
^Kt(^)is s "k ra: *- n tensor history o<t'< + 

X is the spatial point included for any variation 
of other parameters 
"t current time 

By exploiting knowledge of sea ice behavior, (2.1) can be 
reduced to a viscoelastic representation by restricting 
load magnitude and load duration. Further restriction 
results in an elastic formulation. The resultant material 
behavior can be characterized as falling somewhere into 
Figure 2-2. The quantitative limits of the regions have 
not yet been experimentally determined. 

2. 3 Failure 

Category III of figure 2-1 requires knowledge of the 
mode of failure of ice surrounding a circular pile. Most 
investigators work with an assumption of compression failure 
of the ice immediately preceding the pile. However, since 
this failure is quite large with respect to shear or tensile 
failure in simple tests, it is conceivable that tension or 
shear failure may precede compressive failure. 

Additionally, if it is assumed that the greatest level 
of force on the pile will be produced at the point of rupture, 
cracking and large deformation behavior need only be applied 
for postfailure behavior. 

Failure by limiting strain is possible also at levels of 
force below elastic limits and at larger time durations. 
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load 

magnitude 




FIGURE 2-2 

Effect on domain boundaries by load parameters at a 
constant temperature (after Katona and Vaudrey [j8j ) 
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2.4 Material Properties 

There is a large gap in the body of knowledge regarding 
the material properties of ice and nowhere is there a con- 
sistent and universally recognized presentation of values. 
Therefore, for purposes of construction of the model and 
presentation of results of this thesis, non-dimensional 
quantities will be used wherever possible. Some representative 
quantities which may be applied for clarity: 



Quantity 



Typical Value 



Young's modulus 



Poisson's ratio 



Erbp 300,000 psi 

W .33333 



Viscous constant 



Shear modulus 



^ ref x.c-jvjw a xv (psi-sec) 

£>Hp 150 psi 




Elastic limit 



Failure by simple: 



• compression 



• shear 



•tension 



400 psi 
117 psi 
200 psi 



These values are not considered to be rate or time 



dependent for this model formulation, although in general 
they are. 
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CHAPTER 3 

Analytical Modelling 

In the existing literature, two basic mathematical 
models were found for the ice sheet - circular pile prob- 
lem. Both models will be quickly reviewed here and the 
features and limitations discussed. 

3.1 Model of Frederking and Gold 

This model, presented by reference [11], draws its 
basis from the earlier mathematical treatise of Noble and 
Hussein [24]. In this earlier work. Noble and Hussein 
derived the exact solution for dual trigonometric (Fourier) 
series that arise in the solution of a mixed boundary value 
problem, such as that one depicted in Figure 3-1 • 

The assumptions and limitations inherent in the analysis 
of Frederking and Gold are : 

- Interface friction is not considered 

- Plain strain is assumed 

- An infinite elastic medium is used 
Unfortunately, the series presented as the solution is red- 
ucible only for the special case: 

-|; (l- l*') - (l-2w) (3.1) 

(primed quantities refer to properties of the pile) 
Frederking and Gold, however, departed from this analysis early 
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FIGURE 3-1 

Indentation geometry (after Frederking and Gold £| ij ) 
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and used the stress function suitable for the boundary- 



conditions , 






- (l-Z«r»)[r L 



v + 



77]coS©]+D 0 U 2 ["►“*’ a* r**] &*<»**© (3.2 

H«2 



) 



[where is total force on the pile concentrated at 
the center of the rigid pilej 

to derive the radial stress and strain. The infinite series 
of (3.2) was then greatly simplified by assuming a Poisson 
ratio of .5* whereupon the representation for radial strain 
at the interface becomes : 



2Gi t (a.,0) J > cos e - + -5s. 



( 3 . 3 ) 



2-TTcx. * 

The constant D 0 can be uniquely determined in terms of . 

Frederking and Gold then go on to develop a refinement 
of (3.3) in which temperature effects* viscoelasticity* and 
a relation between radial strain and penetration rate are 
discussed. Finally an equation relating total force: 

F*' K(T)LO) 6 0 [“I"] <xK (3 - 4) 

where K(T) is a temperature correction function 

^ is a constant modifying yield stress by 
strain rate effect 

and L(0)is a geometry function based on (3.3) given 

by: 
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here is the "contact" angle 
3.2 Model of Ross, Hanagud, and Sidhu 

Ross, et al . , [27] considered adaptation of one of their 
earlier works [34] concerning elastic-plastic plate impact 
study to the problem of an ice floe surrounding a rigid 
inclusion. The model chosen for this study is based on a 
finite difference scheme for the solution of a plane-stress, 
dynamic problem. The problem is formulated with a two 
dimensional cartesian network. 

The method presented contains the following assumptions, 
although it is clear that the model is not restricted to 
these : 

- pile is rigid 

- loading is by edge dislocation or stress on a finite 
plate 

- perfect adhesion of ice to the pile 

- plane stress situation 

The method of solution of this dynamic problem is 
presented in reference [34], It consists basically of solving 
the continuity equations for the body forces and hence particle 
accelerations. The accelerations are translated into strain 
rates through kinematic relations. Finally the stress-strain 
relations are used to determine stress rates. These stress 
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rates are then integrated numerically into new stress quantities 
whereupon the process is repeated, 

3,3 Criteria for the new model 

For the developement of a model for analysis of the 
circular pile problem, it is felt that the following items 
should be incorporated into the ideal: 

- The analysis diould present variables to include 
stress, strain, their rates, and displacements 

for the complete spatial distribution and at various 
time intervals, 

- The model should at least be a two dimensional 
representation and be capable of plane stress or 
plane strain variation, 

- The model should be flexible enough to accept the 
various material representations of ice including 
non-linear and time dependent ones, 

- Geometrically, the model should be capable of handling 
mixed boundary values, interface friction, and finite 
or (semi) infinite sheet formulation, 

- Fracture or material failure can be easily be incor- 
porated, 

- The model must of course be simple to use and econ- 
omical to apply. 

Unfortunately, these ideals are often mutually exclusive. 

In the next section and in the subsequent chapters, a model 
will be developed that meets at least several of the criteria 
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and will be shown to be useful for application to several 
ice - circular pile interaction problems. 

3.4 Developement of the New Model 

The model to be developed in chapter four is a departure 
from the one described by Ross, et al . , in section 3-2. The 
approach used by Ross seemed interesting in that it is capable 
of meeting several of the criteria required. 

The method of solution remained critically dependent on 
the time integration procedure, however, and any numerical 
inaccuracies or instabilities would propagate throughout the 
solution. The dynamic problem also seemed inappropriate for 
study of slowly varying stress fields. Lastly, Ross gives no 
results for his dynamic problem. 

The model proposed in this thesis, although similar in 
conception to Ross', exploits the radial symmetry associated 
with the physical problem to work in only one dimension. 

The method of solution is greatly dissimilar in that it does 
not treat the dynamic problem and uses a linear system of 
equations to arrive at the final result. The new model 
also treats the case of linear viscous behavior and has 
capabilities of extension into other material behavior. 
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CHAPTER 4 

Developement of a Finite Difference 
Elastic Constitutive Model 



4.1 General 

The equilibrium equations by Wang* reference [31], 
are applied to the element depicted in Figures 4-1 through 
4-4. In two dimensional polar co-ordinates: 



S&r i JL d&re , 6 r ~ 6e _ • 

dr r 5© 1 " r “ 

I c)(3 e i dSjre i 2 / . 

r d© dr + r ^ ^-e 



(4.1) 

(4.2) 



The right hand sides are body forces; the dots indicate 
differentiation with respect to time. 

Since sea ice is a transversely orthotropic material, 
variation of mechanical properties in the Z Co-ordinate 
direction will be due mainly to salinity variation, temper- 
ature profile, and non-uniform thickness. For the purpose 
of this study, however, the ice sheet will be considered 
to be uniform in thickness with zero temperature gradient. 

1 

(Extensions of the present model to account for these effects 
are possible.) For this reason, developement of (4.1) and 
* (4.2) into a two dimensional plane stress model is an 

appropriate choice. 

Before proceding, it is perhaps worthwhile to state 
the assumptions of a plane stress formulation: 
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Stresses acting on a small element ABCD in polar 
coordinates 
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FIGURE 4-2 



Variation of direct and shear stress with angular change 
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FIGURE 4-4 



Sign convention 
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- are zero 'on both faces of the sheet 

- Therefore, Cfej are zero throughout the sheet. 

- The state of stress can be specified by <5r,0e,§Ye 
only and variations of these quantities are 
independent of vertical position. 

4.2 Developement of the Plane Stress Elastic Model 

After Wang [31], the strain-displacement relations 
for small deformations are: 





(4.3) 


** - T + T k U * 






(4.5) 



Stress-strain relationships dictated by the plane stress 
assumption: 



fr = 


T < O r - H6 S ) 


(4.6) 


= 


-g- < &e-^ 6 r ) 


(4.7) 


«re = 


2 ( 1 + w) a 
“ £ °re 


(4.8) 
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Combining (4.1) through (4.8), the following five equations 
are obtained. 



k + 7 h + 


- 0 


(^-9) 


~ do + ^ + ~ 6 re 


-o 


(4.10) 


a T u r 0,.- P0 e ) -O 




(4.11) 


7 u ' + f| u e-| ( 6 9 - 


H6,) = 0 


(4.12) 




2(I + H 0 - o 


(4.13) 



The body forces expressed in (4.1) and (4.2) were 
set to zero since it is not the purpose of this study to 
expose dynamic considerations or impact shock. For the 
range of validity of this model, inertial effects are 
negligible. 

Equations (4.9) through (4.13) (or their cartesian 
counterparts) could be attacked directly in finite difference 
form, but the ensuing number of equations would be quite 
large and difficult to apply. 

If boundary conditions are not to be mixed at a 
continuous boundary, a Fourier simplification is possible: 

Let 6^69 and Uj. be expressed as a Fourier cosine 
(even) series. 

Let be expressed as a Fourier sine (odd) 

series. 



4o 



oo 



<5 r (r,©) - 2Z Or„(r) cos m© 



VvirO 



are similar 



oo 



O^Cr,©) ~ ' O|r0^(i r ) sin in© 



O 



(JIq is similar 

It is then possible to express (4.9): 



(4.14) 



(4.15) 



°o 

L^_, 6 r m Wcos>n©] 

K* = C J 

OO 




i- 00 

(v) sin m©l 

tM=l J 

OO 


6 i-J r) cos^e] 


H") 


^ 0 e (0 cos -m 9~[ 


lM - O J 


i 


— i m ' 



(4.16) 



which reduces to: 



co co 

S& 6 U cos ™ e + C>,. 9 cosmQ 

* mo ^rr 1 ^ 



OO 



+ r ^ j ( ^ " 6 ^Tn ) CO * 



>Vt -O 



(4.17) 



Finally: 



J 



dr ^ r o + 0 + -jT ( 6 



r. - bOf -- 0 
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<* O r| + r Cy 8| + r ( Or, - C>si^ cos© 
[fr^ + r 0 r92 + i(6r,-6 s ^cos2e 



(4.18) 



Thus O^can he expressed as an equation relating only 
like wavenumbers. The set of solutions unique to each 
wavenumber can then be reconstituted into a Fourier series 
by using ah appropriate number of waves. 

Equations (4.10) through (4.13) can similarly be trans- 
formed into an infinite number of sets of five equations: 



A (=) i ]G1 /S. 4. G’r-yy, - 6ewi 

dr + r ° re ™ + ^r 

^7 Ur wi ” ‘g’ ( - V> O e m ) 

( vyi | 

T U r m + ~ “* fT ( 6 evv> “ 7 C>r w ) 



(^.19) 







v* + £r U 



©m ~ lr ^-©wi 



2.(1 + ^ 



>re 



>1 



(Olv^ioo) 
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0© in the fourth row of (4.19) can be explicitly 

expressed as: 

Ce.* = | «ry + u Sh , + w O rw , (4.2o) 

Since no equation contains a derivative of » this 

variable can be eliminated, yielding: 

|;<V(- L ^-)6 r+ + +(^)U 6 

(4- U 6 

( V ‘£ 1 )<5r+ ( O )6r®+C-^— ) +( i Y L )a e 
iu 9 +( O ) 6 r + ('-^)<W(^u r H4-K 

O i m 1 co 

(4.22) 

Then (4.21) becomes the set of linear equations: 

*[i]z.*[a]Z.-Q <*•“ 

where £jj is the 4x4 identity matrix and 



Let the vector r represent: 

^ bu 

6r„ 



z 







(4.21) 
= O 
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where 




expressed by the equation: 



I - P 
Y 




Y 

E 



hi _E - Em 




0 



-_m£ -Em 2 

t 2 - r 2 

wm 

r Y 



(4.24) 



0 



-2(I+f) -m 



(4.23) represents an infinite set of 4 simultaneous first 
order linear differential equations. It can be seen that 
two spatial variables would result in an equation relating 
a derivative in 0 as well. 

The system of equations expressed by (4.23) can only 
be solved numerically. 

4.3 Numerical Solution of the Plane Stress Elastic Model 
The method of finite differences will be exploited to 
solve equation (4.23) 

Considering the finite mesh scheme presented in Figure 
4-5 and in Figure 4-6, the appropriate difference equations 
valid at the nodes are expressed: 

aTlllhi-4Z.-i.HA] Z.-Q 

valid at node 1 (4.25) 



FIGURE 4-5 

Finite difference representation 









1 2 3 4 : 


i 


: • 1-2 1-1 I' 



FIGURE 4-6 



Node numbering scheme 
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tk[i][ 3 Z >i - 4 Z v 1 + Z^ + [a]Z 1 „=o ^.ss) 

valid at the I-th or last node, 
and the central difference 

itfi] { z.- Z-H A 1 1. - 2 <*■«> 

valid at all interior nodes. 

Recasting equations (4.25), (4.27), and (4.26): 



f[Al-i[i]]Z. + l[i]Z.-i[i]Z.=Q 



(4.28) 

(4.29) 



iWZ/WZ.-iWZ.-o 



(4.30) 



Equations (4.28), (4.29)* and (4.30) are then loaded 
into a suitable simultaneous equation - , matrix representation 
as shown in Figure 4-7. 

This representation is then expanded into the full coeff- 
icient matrix shown in Figure 4-8. It should be noted that 
even the simplest scheme, a three node model, will require 
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FIGURE 4-7 

The system of linear equations - before expansion 



47 



1 

2 

3 

4 

5 

6 

7 

8 

9 

10 
11 
12 



'll b 12 b 13 U77 

B 21 B 22 B 23 B 24 ° 

B 31 B 32 B 33 B 34 0 

B 41 B 42 B 43 B 44 0 

— 0 
2h -1 



12 3 4 

" 57 -, B, „ B., , B 



5 

2 



2h 



2h 

o 



o o 

2fi 0 
0 2F 



8 9 

o iL 
2h 

o o 

ooo 
2 

_oo 

h 

0 A 11 A 12 A 13 A 14 ~ 2 h ° 

1 1 



6 

0 

_ 2 _ 

h 

o 



7 

0 

0 

_ 2 _ 

h 



10 11 12 

0 0 0 

-1 



- 1 _ 

2h 

o 

o 



o 

0 

-_ 1 _ 

2h 

o 



A 21 A 22 ^3^4° 2h ° ° 

0 2h 0 

0 ° 2T 



A 31 A 32 A 33 A 34 ° 
A 41 A 42 A 43 A 44 ° 



2TT 0 ° ° ^1 A 12 A 13 A 14fh ° 

-1 J_ 

0 2h 0 0 A 22 A 23 ^4 0 2h 

0 0 2 h 0 
-1 

0 0 0 2 h 

-1 



A 31 A 32 A 33 A 34 ° 



0 

_1_ 

2h 



0 

0 

0 

1 



A 41 A 42 A 43 A 44 ° 



2h 



o o 2h 
A 11 a 12 a 13 a 14 



4(1-1 )-3 


- 


-1 

2h 


0 

-1 

2h 


0 


0 


A n 


A 12 


A 13 


1 

Ai 4 2h 


0 

1 

2F 


0 


0 


4(1-1 )-2 


• 


0 


0 


0 


^1 hz 


A 23 


A 24 


0 


0 

1 


0 


4(I-1)-1 


• 


0 


0 


2h 


0 

-1 


A 31 


A 32 


A 33 


A 34 


0 


0 2h 


0 

1 


4(1-1) 


• 


0 


0 


0 


2h 


A 41 


A 42 


A 43 


A 44 


0 


0 


0 2h 


41-3 


• 


1 

2h 


0 

1 

2h 


0 


0 


2 

h 


0 

2 


0 


0 


C 11 


°12 


°13 


C 14 


41-2 


• 


0 


0 

1 


0 


0 ‘ 


h 


0 

2 


0 


C 21 


C 22 


C 23 


C 24 


41-1 


• 


0 


0 


2h 


0 

1 


0 


0 


h 


0 

2 


C 31 


C 32 


C 33 


C 34 


41 


• 


0 


0 


0 


2h 


0 


0 


0 - 


h 


C 41 


C 42 


C 43 


C 44 



FIGURE 4-8 

Expanded matrix representation 
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reduction of a 144 element square matrix (12 equations in 
12 unknowns). Even a numerical scheme such as the Gauss 
elimination or the Gauss -Jordan reduction scheme is too 
laborious for hand calculation, so a machine solution must 
be employed. 

Figure 4-8 represents the expansion of equations (4.28) 
through (4.30). The elements noted represent the corres- 
ponding position in (4.24). As an example: 



required by forward difference and backward difference 
equations (4.25) and (4.26). The matrix is a square matrix 
distinguished by a band twelve elements wide. Even for a 
small number of node equations the ensuing sparse matrix is 
wasteful in space and computation time. Consequently the 
representation in Figure 4-8 is further compressed to a 
band storage array depicted in Figure 4-9. 

4.4 Non-dimens ionalizat ion 

In order to produce array elements of approximately 
the same order of magnitude, (4,23) was converted to non- 
dimensional form by the following conversion: 




r is evaluated at node 2 



Elements labeled 




represent modification to 




on stress - - — £) 

§REF 




(non-dimens ional ) 
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FIGURE 4-9 

Matrix in the hand, representation 
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RFF 



is taken as the average Young’s modulus through 



the thickness, and length is non-dimensional by 




(non-dimensional ) 



where L = R in or inner radius 

Equation (4.23) becomes non-dimensional and the coefficient 
matrix (4.24) becomes: 



1-P 

Rad 

-vru-* 

RAO 

V>~-\ 



vn - 1 - wi 

RAO RAO 2 RAO 2 

2. - Wi - m 

RAO RAO 2 RAO 2 

Q — 

RAO RAO 



(A-31) 



0 



-2 (In') 



-yy\ 

RAO 



- 1 
RAO 



Here 



RA0= I + 

and will vary within the array according to node position i. 



(4.32) 



4.5 Machine Solution 

Equation (4.23) was implemented in a Fortran language 
subroutine called elastic. For reasons discussed in the next 
section, the number of nodes used was 400. Consequently the 
core storage requirements in double precision amounts to 
approximately 60K words. 
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Solution of (4.23) was by means of an I.M.S.L. 
subroutine called LEQT1B [16]. This method of solution 
of a system of linear equation is described by [23] as 
the Crout factorization. The Crout factorization is similar 
to Gauss reduction methods but the lower triangular decom- 
position is used instead. Partial pivoting and row normal- 
ization are employed. 

LEQT1B was slightly modified into a double precision 
routine and labeled dpband by the author. This routine is 
included as an integral part of subroutine elastic and is 
listed in Appendix A. 

Subroutine elastic was written for the multics system 
implimented on a Honeywell 6l80 used by -the Massachusetts 
Institute of Technology Information Processing Center. 
Approximate CPU time for one solution was 1 6 seconds. 

A flow chart for subroutine elastic is presented in 
Figure 4-10. A symbol table, description of common blocks, 
listing, and a user's guide are included in Appendices A-l 
through A -4. 

4-6 Numerical Stability 

The number of stations chosen is 400. Table 4-1 and 
Figure 4-11 illustrate the convergence of the numerical 
solution as the number of stations is increased for the 
particular case of zero wavenumber. Having checked the 
convergence of the numerical solution, the question now is 
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Subroutine Elastic 




Set material constants 



Zero arrays A, B 
Zero arrays C, D, R 



Derive fourier coefficients 
for boundary conditions 



Display fourier coefficients 
of boundary conditions 



ij 



Read wave numbers desired 



Load coefficient matrix A 
Load boundary conditions into 
coefficient and solution matrices 
Solve 

Store wave solutions 

Calculate stresses, strains, and 

present output 

Combine series 

Zero arrays A, B 



I Q2/t\ 



Zeroabx 



Zerocdr 



Loader 



Courier 



WritrO 



Loada 



Loadbdy 



Dpband 



Strain 



Recomb 



Zeroabx 



wavenumber = wavenumber desired? 



Display total solution 




Writr2 









Stop 



FIGURE 4-10 

Block diagram for subroutine Elastic 



Simps 
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Table 4-1 



Wavenumber 



Maximum Oscillation 
(per cent) 



0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 



0.0075 

1.7500 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 



Radius - outer = 20. 

Radius — inner = 1. 

u^ = Uq = 0. inner boundary 
( 3 p = -.66666667 outer boundary 
6 f 0 = 0 * outer boundary 



Maximum oscillation in numerical solution by wavenumber 
(solution carried to five digits) 



5b 




Effect of interval spacing on numerical stability - radial stress 
(effective maximum swing - representative points plotted) 
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whether the numerical scheme is consistent, i.e. whether 
this convergent numerical solution is the true solution. 
This is done in the following section. 

4.7 Comparison with Analytic Solution and other Tests 

Fortunately, the analytic solutions to the relations 
expressed in equations (4.21) can be obtained and are 
presented in Appendix B. The analytic solutions and their 
computed counterparts are illustrated in Figures 4-12 
through 4-l4 for wave numbers 0, 1, and 2 respectively 
for the case of radial stress at R out/R in = 20. 

Test A This was run to determine the effect of radius 
on the inner stresses. It would be expected that for a 
stress loading of the outer boundary, the importance of 
the radius on the inner boundary would decrease with 
increased plate diameter. This has been determined to be 
so numerically, and some results are presented in Figures 
4-15, 4-l6, 4-17 » and 4-l8. It is also noteworthy that the 
higher wave numbers have very little effect on the inner 
stress solution. A conclusion here can be drawn that will 
save much computational effort if only a limited amount of 
data is desired. 

Test B This was run to determine the sensitivity 
of the solution to the Poisson effect. The results for 
four wave numbers are presented as Figure 4-19. Wave 
number one appeared to be almost entirely insensitive, 
whereas wave number zero appeared to be the most sensitive. 
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stress 




Numerical and analytical solutions wavenumber zero 
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stress 
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stress 







Effect of radial size on radial stress distribution - wavenumber 
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FIGURE 4-15 
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-3.0 



STRESS 



-2 

x10 



6 r =-666667 
ROUT = 20 

5 

WAVENUMBER 




-5.0 



0 



1 




Variation of radial stress at the inner boundary 

as a function of poisson ratio for several wavenumbers 
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The next two tests were of load magnitude and shape 
and designed to test for linearity. 

Test C The load magnitude was varied by integral 
multiples, retaining the same distribution. Results were 
exact multiples as well. No results are presented. 

Test D The test of load shape was more of a test 
of the Fourier subroutine which converts boundary condi- 
tions into series. The following radial loading was 
supplied : 



Or 

U r 



£L 


#2 


#3 


a 1. cos 0 


= 1.-1. cosO 


=1. 


i - o. 


=o. o 

= 0. ? 


=0. 


INKER ~ O. 


=0. 


INNER - O. 


= 0. = 


--0. 



The results were excellent and are presented in Table 4-II. 

Test E The last test concerned the radius at which 
the finite plate could be considered infinite for practical 
purposes. This test is similar to Test A. This test was 
somewhat inconclusive since above a certain radius some 
solutions became unstable again and innacurate. A practical 
limit of radial spacing was chosen at n=400, since the 
marginal accuracy of finer differences was outweighed by 
increases in computational time, cost, and round off error. 
Results of this test are presented as Test E in Table 4-III. 
4.8 Comparison with Ross* Solution 



In Section 3*2 a review of the analysis of Ross was 
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Table 4-II 



Radius 


Run 1 


Run 2 


Sum 


Run 3 


1.00 


9.72992 


-8.2323 


1.49762 


1.49759 


2.90 


5.34873 


-4.2905 


1.05823 


I .05824 


4.80 


3.36513 


-2.3445 


1.02063 


1.02047 


6.70 


2.46372 


-1.4538 


1.00992 


I.OO 989 


8.60 


1.95948 


-0.95387 


1.00561 


I.OO 551 


10.50 


1.64134 


- 0.63805 


1.00329 


1.00328 


12.40 


1.42509 


-O .42309 


1.00200 


1.00200 


14.30 


1.27047 


-O .26917 


1.00130 


1.00119 


16.20 


1.15562 


-O .15491 


I.OOO 65 


1.00065 


18.10 


1.06816 


-0.067790 


1.00037 


1.00027 


20.00 


1.00002 


- 0 . 11921 e -04 


1.00000 


1.00000 




Up j Uq = 0 

inner 


u rj u e =0 

inner 




u rj u© = 0 

inner 




outer 

6 r =COS0 

6©=0 
Rout- 20 


outer 

6 r = 1-cosG 

6©= 0 
Rout=20 




outer 

6 r =1 

6©=0 
Rout = 20 



Test of elastic routine by variation of load 
shape 
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Table 4 -III 





outer 10 
inner 


5 % deviation Instability 


Wavenumber 


©Bout 


@Rout 


@Rout 


0 


>100 


>100 


15 


1 


>100 


50 


10 


2 


>100 


> 100 


75 


3 


25 


>100 


>100 


>4 


<25 


>100 


>100 



Practical extent of finite condition at certain wavenumbers 
by variation of outer radius 



67 



presented. This paper, reference [27], also included an 
example in which the loading depicted in Figure 4-20 was 
analyzed. Below the stated elastic limit, it was possible 
to compare the model developed in this chapter. An effort 
was made to correlate this result with a similar loading 
applied to the method developed in this chapter. The 
comparison is presented as Figures 4-21, 4-22, 4-23. 

Although the basic shape and critical values of cross- 
over are similar, discrepancies were noted with respect 
to the hoop stress and also with respect to the order of 
magnitude. By changing the displacement loading magnitude 
to .01 inches vice .01 feet, the numbers were much alike. 

No explanation was found for the difference in hoop stress, 
but from equation (4.20): 

<Jr M + ^ U 6m + 

Evaluated at • U 6) U r = 0 (boundary conditions) 

produces 

or a non-zero value for any non-zero radial stress. 

This loading was also run with 3* 5* and 11 wave numbers 
and no effect was found on the value of inner face stress, 
although significant variations existed toward the edge. 

This is as expected, since the influence of the trigonometric 
components wi)| do not substantially contribute to the center 
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FIGURE 4-20 

Discontinuous edge loading of a finite plate about 
a rigid inclusion (after Ross, et al , [27]) 
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Radial stress vs. angular position - comparison vfith Ross' solution 
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Hoop stress vs. angular position - comparison with Hoss* solution 
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Shear stress vs. angular position - comparison with Ross' solution 
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stress as the edge of the plate is extended toward (perceived) 
infinity. 

Figure 4-24 displays the radial loading represented 
by three Fourier coefficients and the original loading. It 
can be seen that even for this discontinuous loading, three 
wave numbers presents a good approximation. 

4.9 Peculiarities and Range of Validity 

This elastic formulation is derived from the basic 
plane stress equations and consequently should be useful 
where plane stress is appropriate. The model achieves 
significant accuracy up to a ratio of outer radius to inner 
radius of 25:1 • Beyond this, the solution diverges from 
the analytical derivation. 

The worst case is at wave number one (Figure 4-25). 

It is doubtful that the infinite case can be modeled by a 
finite difference technique in this problem, since a 
boundary loading for the first wave number at all radii 
will produce a significant solution at the pile. Consequently, 
this elastic model is useful for radial limits of outer radius/ 
inner radius ratios less than 25. 

The buckling limit of a thin plate is an important 
consideration. The analytic derivation of one buckling load 
is presented as Appendix C. The equations of plane stress 
do not allow for deflections in the third dimension, therefore 
they will not predict the failure by this mode. External means 



ifWOOCMr-O lO CO CM 
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FIGURE 4-24 

Displacement loading (radial) applied to boundary 
of Ross' solution - comparison with first three 
wavenumbers 
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must be applied to the program to preclude this occurance. 

The plate is assumed to be attached to the rigid in- 
clusion (frozen-in problem). The partial movement of either 
boundary cannot be supported by this method, although it is 
conceivable that a more general two-dimensional model could. 

Finally, the model can support an application of forces 
or displacement at either boundary, which can in turn be 
discontinuous. The degree of refinement is limited by the 
wavenumbers considered. The accuracy is limited by the 
circumf erential spacing chosen. 
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CHAPTER 5 

Viscous Extension of the Model 



5.1 Genera l 

Considering the simple rheological model in Figure 5-1 > 
the following differential relations hold: 



• • 

£r ” £r viscoos 


+ 


• 

( r ELASTIC. 




(5.1) 


• • 

6© — 6^ Vis coos 


+ 


« 

£ e ELASTIC 




(5.2) 


# • 

(p© = 6p©viscoos 


+ 


• 

£ re elastic 




(5.3) 


• 

where £ 


V1*C003 - ^UNCTION 


(6) 


(5.4) 


Maxwell model is 


chosen for its 


simplicity 


and 



because it is possible to generalize into a non-linear 
creep law. Additionally, some investigators have modeled 
sea ice behavior this way experimentally, and physical 
values are available [33]. 

5.2 Incorporation of the Viscosity Relations 

The stress-strain relations of (4.6), (4.7 ) » and (4.8) 
can be expressed using (5.1), (5.2), and (5.3): 



( 0 r - O e ) + -g- ( - v^6 e ) (5.5) 

(© » ( O© - W Or ) + -gr ( 6 e - ^ 6 r ) (5.6) 



(5.7) 
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Model for Maxell liquid (after [j3] ) 
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These relations can he incorporated into a series of equa- 
tions as in section 4.2 yielding: 

<5 - 81 

Where : o £ £ oo 

[A] is the same matrix as expressed by (4.24) 

And 






o 


1 

0 


2 

- E Wi oi 




r 1 


r 1 




-E\w* 

° 


c 2 1 


o 


— fc Vv\ ot 


o 


Q E«(' + ^ 


Evwoc(h-»>) 






r 


o 


2 01 0+w) 0 


O 



(5.9) 



5.3 Solution 

Numerical solution of the series of differential 
equations (5.8) is carried as before, only the unknown 
variable, ~7 , represents rate at time "t • The value 

Z — * 

used on the right hand side of (5.8) is the elastic 
solution for y' . Note that boundary conditions of load 
rate must be employed. The solution yields: 



A. 




to 



= r 



( 5 . 10 ) 
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This differential equation in time can he solved using the 
Euler or "one step" method to arrive at a solution for the 
next time step, where the process is repeated. 



. t t *fc + A-fc 

(ZJa* + z„= z* 



( 5 . 11 ) 



Non-dimens ionalizat ion with respect to another variable 

is required. The variable chosed was oP , the viscosity 

-I 

coefficient, with units of LTM- 

If ©Preference is chosen, and mass is preset by consideration 
of stress non-dimensionalization (Section 4 . 4 ), then derived 
time : 




E 



REF RtsF 



( 5 . 12 ) 



Therefore any value of non-dimensional time is related: 




TpHYSlCAL 

”T*REF 



( 5 . 13 ) 



As an example, using typical values for sea ice [ ], 

El ref = 2083 p'sf 

^ref = 1.2788* lO' 1 (psf-set)' 1 

implies 

E = i. 



and 



cP 



Thus 



Tref - 2607 sec 



or one unit of non-dimensional time 
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A flow chart of the numerical computer solution and 
time integration is presented as Figure 5-2. 

5.4 Convergence in Time 

The Euler method chosen for time integration assures 
convergence if the time step is small enough. Since it 
is not possible to analytically determine the solution, as 
was done for the elastic case, one approach is to run a 
general solution for increasingly small time steps until 
no numerical improvement is achieved. Unfortunately the 
solution time can be prohibitive. The results of two time 
convergence test series, for wave number zero and one, were 
found to converge at a time step of .001 T , beyond which 
significant improvement was not achieved. The time increment 
to achieve this numerical stability comes out to be a constant 
independent of stress, for the upper bound prediction. 

5.5 Range of Validity 

The viscoelastic formulation presented in this chapter 
is not valid where the elastic theory is inappropriate or 
where the numerical limitations of subroutine elastic do 
not apply. The range of validity is further restricted 
on how well and on which range of parameters ice can be 

1 

modeled as a linear viscoelastic material. 
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Subroutine Visco 







FIGURE 5-2 

Block diagram for subroutine Visco 
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CHAPTER 6 

Analysis of Ice - Circular Pile Behavior with the Model 



The previous chapters have dealt with the need for an 
analytical tool and the construction of a mathematical 
model* Without exhaustively analyzing the many combinations 
of loadings and important parameters that affect the circular 
pile problem, some interesting results will now be developed. 
6.1 Material Behavior 

Katona and Vaudrey [18] have presented a logical seq- 
uence of obtaining ice material behavior experimentally for 
a particular geometric configuration. The results of these 
data are then fitted into curves which can then yield a range 
of validity of ice material law as discussed in Section 2.3 
and with Figure 2-2. 

Of course, since this model formulated in this thesis 
responds only viscoelastically, the transition from visco- 
elastic to viscoplastic behavior cannot be obtained. The 
model can, however, define the elastic -viscoelastic transition. 

As an example, consider the loading situation depicted 
in Figure 6-1. Physically, this could represent the rigid, 
frozen-in pile loaded by a large, thermally expanding ice 
sheet. Neglecting any temperature effects to viscosity and 
elastic modulus, the loading curves of Figure 6-2 are con- 
structed. Stress-strain diagrams similar to Figure 6-3 show 
the viscoelastic creep behavior. The apparent relaxion of 
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Viscosity coefficient oi = 1 . 

Elastic Modulus E = 1. 

Non-dimensional quantities, 
see text in Sections 4-4 and 5-3 




FIGURE 6-1 

Circular ice sheet - rigid pile geometry and loading 
configuration used in viscoelastic "behavior analysis 

of Section 6.1. 



m 




FIGURE 6-2 

Viscoelastic behavior - Effective strain vs. load duration 
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Effective , 

Stress xlO - 1 



&ef f = ( ^ b r + b e“ b r b © + 3b r © ))2 




Viscoelastic behavior - stress-strain diagram showing 
creep for loading condition of Figure 6-1, at the pile- 
ice sheet interface. 
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elastic modulus is then limited to a value of 10$ decrease. 
This then becomes the transition boundary between visco- 
elastic and elastic behavior of Figure 2-2. 

In the example, loadings of those shown in Figure 6-2 
are used to define a region of validity in Figure 6-4. 

Using actual values such as those found in Section 
2-4, it can be seen that a time of approximately 10 minutes 
defines the region at low stress levels (below yield). 

The effective stress and strain plotted in Figure 6-3 
are defined : 



For this loading condition, it can be seen that a 
value for modulus of elasticity is about one half the uni- 
axial value. This will vary as well along each point in 
the radius of the ice sheet. Interestingly enough, the 
relaxation to 90 $ of its former value occurred at the same 
time, despite the variation. Under a more general loading, 
not confined to the zero wavenumber, it is likely that the 
material will behave not in the same region of Figure 2-2 
over its entire geometry. 

6.2 Under Rate Loading 

The example in the previous Section was subjected to 
rate conditions on the loading. This was done to test the 
capacity for this parameter in the model. The load-time 
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Load t 

Magnitude xlO 




FIGURE 6-4 



Load duration - non-dimensional time 



Viscoelastic behavior - Elastic-Viscoelastic transition 
from' load curves of Figure 6-2. By definition, relaxation 
of elastic modulus is limited to 1C$. 
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curves for the loading of Figure 6-1 are presented in Figure 
6-5 for three loading rates. The greatest rate is very 
close to the maximum of 5000 psi/min suggested by Katona 
and Vaudrey to avoid dynamic effects. No general conclusions 
are drawn here since the size of the initial load and the 
various rates have little significance without a specific 
physical problem to simulate. 

6.3 Analysis of a Friction Coefficient 

The loading situation in Figure 6-6 was applied to the 
model to determine the ratio of shear stress to (compressive) 
radial stress at the rigid pile interface. This could be 
construed to be an indication of a static type of frictional 
coefficient, but since the model also assumes capability 
of the interface to support a tension, in this region "fric- 
tion coefficient" would be meaningless. 

The analysis proceeded with those values that this edge 
loading produced. Defining the coefficient of static friction 
as : 




(6.3) 



it can be seen that this coefficient will vary about the 
interface at least to the point at ® s ”j5T w ^ 9rs becomes 
tensile and a coefficient of friction becomes meaningless 
(Figure 6-7). 

At a value of © * .31 radian (given this loading). 



= .3 
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Visco-elastic behavior - vs. time curve under 

three loading rates, at the pile - ice sheet interface. 
Loading condition of Figure 6-1 applies, with 6p=. 1111x10“^ 
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Rout =20. 

2 

£>P =-.333 cos © outer boundary 

b|- 0=-.333 sin© cos© outer boundary 
for 

and 0 . otherwise 




at the inner boundary or pile-ice 
sheet interface: 



Rin =1. 

Up =Uq = 0 . considering "friction" 
and 

bp @ =u r = 0 . "frictionless " 



FIGURE 6-6 

Loading condition circular pile - ice sheet interaction 
for cases discussed in Sections 6-3 and 6-4* 
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Ratio 




FIGURE 6-7 

Ratio of shear stress to radial stress at the circular 
pile - ice sheet interface, loading of Figure 6-6 applies. 
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Now if there is no ice - pile adhesion present 
and if the frictional coefficient of smooth sea ice on 
a smooth steel circular pile is about .3 as some inves- 
tigators indicate, then it is clear that for regions 
beyond .31 radian CC© will be nonzero. A situation in- 
volving mixed boundaries and separation will than be evident. 
By changing the loading to: 

“ Oy. Co^ © 

Or© - - svv\ 0 

it was found that the coefficient of friction at the inter- 
face changed only slightly. Since for the wavenumber zero, 
only normal stresses are conveyed (being an odd function) 
and the higher wavenumbers lose their influence with increased 
radius, frictional forces are important only in the region 
|0| £ .31 radian for non-axisymmetric loads. 

A conclusion can be drawn that for the model of the 
circular pile - ice sheet interaction and unless ice adhesion 
is present, frictional behavior is confined to a small area 
preceding the pile and may be unimportant in analysis. 

6.4 Analysis of the Behavior of Ice Floe - Pile Interaction 
Comparing Adhesive and Non-adhesive Behavior 
Since the influence of friction at the circular pile is 
minimal, it is necessary now to compare the behavior at the 
interface using complete adhesion (frozen-in situation) and 
and using complete circumferential freedom. 
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The loading of Figure 6-6 again was applied to the 
model to consider the effect on the radial stresses of a 
"frictionless" boundary. Since the loading is the same 
for both, the net force felt by the pile is the same in 
both cases. Figure 6-8 presents the variation of radial 
stresses in both cases. The point where adhesion is 
required is precessed to 1.48 radians. 

Noble and Hussein [ 24 ] obtain a value of 1.24 radians 
analytically in their formulation. Severe differences do 
exist in the problem as their model does not support tension. 
Therefore it would be expected that in the new model, 
the radial stress on the compression side would be supporting 
the lost force on the tension side, with basically the same 
shape in the loading curve. 

At any rate, a conclusion can be drawn that the friction- 
less case is the worst as far as local normal pressures on 
the pile structure are concerned, and that the maximum force 
can be as high as three or four times the frozen-in case. 

Of course, the total lateral load remains the same. 

6.5 Locations of Maximum Stress and Strain 

As a final use of the model, information from the loading of 
Figure 6-6 will be presented. Physically, this edge loading 
might be produced by a wave driving force impinging on the 
circular floe. No failure mechanism is included, since 
some investigators feel that failure is ductile, caused by a 
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Stress normal 
to the pile xlO - ^ 




Comparison of normal (radial) stress at pile - ice sheet 
interface, with and without friction. Loading of Figure 
6-6 applies. 
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limiting strain, and no reliable theory is available. 

Therefore no estimation of a maximum total force will be 
offered, but a good representation of stress and strain 
distribution can be presented. 

Figure 6-9 depicts the variation of normal (radial) 
and shear stresses throughout the ice sheet and at selected 
angles © . The greatest radial stress is at 0=0° as 
would be expected. The crossover to tension occurs at 0=90°. 

The maximum shear stress occurs at about 5^° again at the 
inner face. 

Figures 6-10 and 6-11 show the variation in vertical and 
hoop strain around the pile. If tension strain failure is 
the failure mechanism, it can be seen that failure will occur 
in the sheet vertically (a horizontal crack) and horizontally 
(vertical crack) at angles of 0° and 5^° respectively. Shear 
strain follows the shear stress and in this case is 2/3 the 
§r© value. Shear strain is shown in Figure 6-12. The values 
of these strains are plotted in Figures 6-10, 6-11, 6-12 only 
at the interface, since this is the location of their maxima. 

It should be pointed out that the basic shape is maintained 
some way into the ice sheet, so that a crack x*ould be init- 
iated at the ice - pile interface and propagate radially 
outward. 

It would seem heuristically that the cracks at 5^° 
would occur first due to the greater magnitude of the strain. 

It is then possible that the crack at 0° will then follow due 

the redistribution of stress. Both types have been reported [15]. 
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Elastic behavior - Loading of Figure 6-6 - Variation 
of normal (radial) and shear stresses throughout the 
ice plate 
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Elastic behavior - maximum vertical strain location on ice - pile 
interface. Plot of vertical strain vs. angular position at radius 
equal to 1. Loading of Figure 6-6 applies. 
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FIGURE 6-11 

Elastic behavior - maximum hoop strain location on ice - pile 
interface, plot of hoop strain vs. angular position at radius 
equal to 1. Loading of Figure 6-6 applies. 
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Angular position in radians 



FIGURE 6-12 

Elastic behavior - maximum shear strain location on ice - pile 
interface. Plot of shear strain vs. angular position at radius 
equal to 1. Loading of Figure 6-6 applies. 



100 



CHAPTER 7 

Conclusions and Recommendations 



7.1 Conclusions 

The model developed in the preceeding chapters is 
applicable to a circular elastic or viscoelastic sheet 
surrounding a rigid circular inclusion. The conditions of 
plane stress are assumed, although plane strain could he 
very easily implemented. Loading can he hy displacement 
or stress at either boundary, although mixed conditions 
at a particular boundary are not allowed. The model is 
favorably disposed for inclusion of failure criteria, but 
none are included or proposed. Tension is required at the 
interface to support a valid operation and this corresponds 
to a frozen-in situation. This model can accept time dep- 
endent loadings below dynamic significance. 

The model is very thrifty with computational resources 
and produces a very high degree of correlation with existing 
analytical solutions. The model is capable of considerable 
expansion. 

Some results with a fictional loading predict frictional 
effects are present only at a small region immediately pre- 
ceding the pile. If angular movement is allowed, the resul- 
tant loss of shear requires an increase in normal stresses. 

If a separation of boundary is to occur, and occurs where 
tension cannot be supported, then the contact angle will be 
1 90°. 
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The range of time where ice can be considered elastic 
(below yield levels ) with a uniaxial load is : 

.. c .2$ 

time < 

Other loadings could be analyzed. Rate of loading is 
important insofar as time to failure and required force. 
Finally, viscoelastic behavior predicts ice failure (and 
maximum forces) for loading well below immediate fracture. 

Lastly, it is possible to qualitatively see failure 
by limiting strain or shear at locations in the ice sheet 
that are physically relevant. Without knowing the values 
of biaxial strengths of the ice material, it is not possible 
to predict a maximum force before fracture, but it is likely 
that a limiting strain or shear strength precedes ultimate 
crush failure. 

7.2 Recommendations for Future Work 

The model developed in this thesis is extremely rudi- 
mentary and can be used for pile - ice sheet interaction 
under only the most generous conditions. Improvements can 
be made, however, to increase the usefulness into problems 
with more physical application. These improvements by order 
of importance are : 

(1) Driving Forces 

Ice sheets and ice floes are driven by the effects of wind 
and current shears. It would be desireable to develope an 
analytical or numerical procedure whereby these shears can 



102 



be approximated into edge loading conditions. The size of 
the ice sheet considered in this model would require negli- 
gible wind or current shear itself, so that the two-dimensional 
approximation could remain valid. 

(2) Failure 

The failure mechanism or mechanisms can be easily applied to 
the model if a good choice could be found. To date, inves- 
tigators have proposed Tresca or Von Mises type criteria, 
failure by limiting strain, rate dependent crush strength, 
among others. The behavior of ice is poorly understood in 
biaxial loading and experimental work needs to be carried 
out to properly extend the model into realistic and workable 
significance. 

(3) Two-dimensional Extension 

The model should be given another co-ordinate in the plane 
such that more general geometries can be analyzed. Of course 
this could mean larger quantities of data to be handled, and 
an attendant loss of simplicity. The problem of the mixed 
boundary conditions and different shapes of the ice floes or 
structures could be adjusted to fit physical requirements 
with a high degree of confidence. 

(4) Multiple Structure 

The model, once in a more general two-dimensional state, 
could be extended to study the effect of multiple structures 
such as support pilings on a pier or multilegged towers. 
Possible reductions in forces by interference could be shown. 
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(5) Other Material Behavior 

Once in a two-dimensional form, the model could be extended 
into non-linear and non-elastic behavior in a straightforward 
manner. The one -dimensional form presented in this thesis 
is valid only where superposition is valid. 

(6) Other Structural Behavior 

Experimental work with the model could be extended into 
non-rigid analysis for study of structural responses and 
forced vibrations. 

Data on the material properties of ice in a biaxial 
state of stress need to be acquired for realistic imple- 
mentation in the model. Once this and the six stated improve- 
ments of the model are implemented, the behavior of the ice 
sheet-structure problem can be fully analyzed and more general 
design equations or curves offered. Additionally, the 
influence of other important parameters such as temperature 
profile or variation of ice thickness can then be included 
or discarded. 
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APPENDIX A-l 
Program Variables 


Variables in Common 


Purpose 


n 


number of radial spaces 


nn 


number of circumferential spaces 


nnpls 


nn+1 


ind 


indicator of boundary conditions 
used (see listing under load- 
bdy) 


m 


wave number 


ru 


poisson ratio 


el 


Young's modulus 


Pi 


3.1415... 


alpha 


viscous coefficient 


mdes 


number of wave numbers to be 
processed 


rout 


outer radius 


rin 


inner radius 


idsp 


4xN t number of radial display 
intervals 


ibar 


termination variable 



Variables Consistently Used 

irowmax number of simultaneous equations 



to be solved 
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tstop 


stop time limit (visco elastic 
solution) 


tstep 


time step interval 


time 


time since start (elastic: time = 0) 


disint 


time display interval of viscoelastic 
solutions 


streslim 


elastic limit of stress 


strnlim 


limiting strain 


rad 


radius 


hoop 


hoop stress 
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APPENDIX A-2 
Description of Arrays 

Array Subroutine 

/blokl/a^ l 6 oi 3 17 ) elastic , zerabx, loada, strain , 

double precision loadbdy , visco, dpband 

purpose: holds coefficient matrix in band compressed 
form - after call to dpband holds the factored matrix 
/blok2/b(l601,l ) elastic , zerabx, loadbdy , 

visco , loade , boundb , strain , 
dpband 

purpose: solution vector of the system of linear 
equations represented by a 

/blok3/£ (204,21) elastic, zercdr , recomb, 

writr2 

purpose: stores total recombined solutions at every 
eighth node (row) and all circumferential positions 

/blok4/d(ll,8) elastic , zercdr , fourier , 

loadbdy , writrO 

purpose: stores fourier coefficients of boundary 
conditions for 0 through 10 wave numbers (row), inner 
conditions (first four columns), outer conditions 
(second four columns) as determined by variable ind 
rl(8,21) loader , fourier 

purpose: stores boundary information from main before 

fourier decomposition 
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ws(21) fourier, simps . 

purpose: working vector to be integrated by simps 
/blok6/xl (14436^ elastic , zerabx , dpband , visco 

double precision 

purpose: working storage used by dpband 
/blok8,£r(204,ll) elastic , zercdr , recomb 

purpose: stores elastic wave number solutions for every 
eighth node 

/blok9/t (51*11) elastic , zercdr , recomb 

purpose: stores hoop stress wave number solutions for 
every eighth node (similar to /blok8/) 

/bloklO/tt (51,21) elastic , zercdr , recomb , writr2 

purpose: stores total recombined elastic hoop stress 
solution (similar to /blok3/) 

/blokll/rr(l6o4,l) visco , loade 

double precision 



purpose: working vector 
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APPENDIX A -4 
User Notes 



Program Operation 

The program flow is described in Figures 4-10 and 5-2. 
The material constants E, v > , ©t are set in either the main 
calling program or subroutine elastic prior to compilation. 

In the listing presented in A-3, E and o< are non-dimensional 
with a value of 1.0. The program is not restricted to 
non-dimensional solution, however, and any consistent values 
could be used. 

User inputs are via keyboard and consist of 3 inputs. 

For elastic: 

1) mdes - after reviewing fourier coefficients of 
boundary conditions, the user is asked for the 
number of fourier coefficients to be used in the 
solution. This feature is incorporated to save 
calculation of meaningless higher wave number 
solutions. 

2) The user is queried on whether a viscoelastic 
solution is required. This feature is used to 
permit an elastic solution only. 

3) Display of theta information - This feature 
automatically combines elastic wave solutions and 
calculates circumferential variation for intervals 
of .05 TT and radial variation of 50 intervals. 
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Display here is set for .01 7T and radial variation 
o£ ten intervals. 

If a viscoelastic solution is desired, additional 
keyboard inputs are : 

1) tstop - desired stop time, entered in units of 
non-dimensional time. 

2) tstep - interval calculation time - This is 
important as a small step must be chosen to ensure 
convergence. 

3) disint - time display interval of wave number 
solutions - It is anticipated that not all small 
discrete time intervals will require display. 
disint should be an integral multiple of tstep . 

Loading Input - rout must be set in the main program. The 

vectors putl and put4 contain boundary condition loading 
information as determined by ind. also lent from main. 
putl - put4 contain diecrete data for loading 
circumferentially through © = TT . (The other half 
is identical by symmetry). 

Inner Boundary Conditions Outer Boundary Conditions 
ind putl put 2 put 3 put4 



1 u*. u© 

2 6 |* 

3 U. r U© 

4 6y- 6r© 

Caution must be exercised with ind 



U.y U© 

6 . 6 ,. 

U*- U.9 

= 2 y since a rigid body 
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translation may be produced. 

Output of Program - Three distinct listings are produced: 

1) Fourier coefficients of boundary conditions - all 
wave numbers 

2) Wave number solutions of , (5 r © ,U. r , U.q, 0© , 0 CFF 

3 node, radius for upto 400 nodes. In the 
listing of Appendix A-3, only ten radial intervals 
are displayed, idsp will vary display points. 

Wave number solutions are displayed for the elastic 
case and at each time display intervals determined 
by disint . The headings "stress" and "strain" 
contain an effective stress and strain as defined 

by *‘ I 2 - 

6 ( 6 r + C>£ " + 3C>re) 

€ ~ ( ^ ^0 * 

3) Stresses by angular position - contains up to 
fifty radial intervals with the total recombined 
fourier series representation of the elastic 
solution displayed at intervals of .OlTT . 

Other Notes 

Viscoelastic boundary conditions (rates) are preset 
to zero in this program, but are not necessarily restricted 
to this value or even a constant value. Rate information 
must be in size with integration of time step and convergence 
must be ensured overall. 
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streslim and strainlim are preset constants corresponding 
to an elastic limit and a strain failure limit. The solution 
is deemed inappropriate beyond these levels and a solution 
is terminated for that particular wave number. 

The accuracy of the fourier subroutine is dependent on 
the circumferential interval. For complex loading it is 
conceivable that this interval can be decreased by increasing 
nn and of course pertinent array sizes. 

Storage of solutions in the B, C, and R arrays for a 
particular variable are in every fourth row, e.g.; 

Z.*b(cj),b(c + lj),b(u2j) ) b(u3j) 



solution 
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APPENDIX B 

Analytic Elastic Solution for the Case of a Circular Ring 
Wavenumber m=0 

The stress function is known to satisfy the biharmonic relation 



o 



(B.l) 



In a two-dimensional elastic material. For the case of no 
circumferential variation (as in wavenumber 0), (B.l) reduces 



to 


(£ + !i_Y^ + l *±) r0 


or 


aH i_ dH _ j_ oH , a't* 

Tr"* + V dr 3 r l dr 3, * V 3 dv- 


Let 


\* 


Then 


d^f _ it __ 1 at 

el *■ tl^elv v* <1 ^ 




l± u A » /dH <m 

dr 1 "3r\ e\v-/ V 1 V d d \ ) 

aH 1 / elH _ d l «f 2 a 4 ** \ 
2 dV* Z lL\) 

. i /aH / aH , ,, aH _ 


(B.2) becomes 


dH. 4 dH + 4 al± _ 

Tp 4 d r 


The general solution of which is 




* » e*' + C z e^ + !f C. 3 v ' 



*t) 



(B.2) 



i- - c, (\©^ e* , '° v + + C-4 






or 
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Finally * * C.^ 1 lo^v- + C z lo^v + C 2 v- 1 * C 3 * C 4 

H°w 6 f = V 3P * C, (U2l» v ) +1C» ♦ ^ (B.3) 

6 e = -0 - C-. U + iVWc, - ^ ( B - 4 ) 
0 fe * o 

Strain relations 






do 
d v* 



In plane stress 



6* -- o 




*r® 



dv _ 

dv - V- 



do 

dv 





(B.5) 




(9 



- H 




(B.S) 



Integrating (B.5) 

Eu * ^(B.3)^ - »* $(B.4)dr 

Into Eu. - C, [* ('-‘5**) * Z(l-*X* U^v-v-)] + 2(^(1-^ 

-C 3 (»+^~ + C 5 (B.7) 

(B. 6 ) becomes 

E a -- (3-w) +2v(l-‘ J Mo< i y] v 2C*(l-p)v* 

- C 3 (up') 7 (B. 8 ) 

For (B.7) to have the same representation as (B. 8 ) 

C v - C* -- O 

Therefore, for wavenumber 0 

Eu. * 2 Ci (\-p)r- C-jU+w^-Jr 

6,. - ^ c.i * ^ 

6. • 2 C t - ii 



(B.9) 

(B.10) 



Wavenumber m=l 
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Prom the general solution for the stress function in circular 
coordinates presented by Timoshenko [29], that portion cor- 
responding to the first wavenumber is : 

*f* a t © so* $ + ( y 3 + a', r ' *■ fe, vIc>y ) cos 9 



Since 



- J_ if 1 

1 5 © l 



<3.> tl 
3v- 



Then 



r - 3 /J_a4\ 

' zZ u iej 

6*. * [laA)v 1 + 2b,** -2 c*.,*- 3 ] cos © 
6© - (6^v- 4-2a.',v- 3 +- V>) v’ l )cos© 
6 rfe ' ( 2 b, ir - 2 A.| *" 3 + \>\ V ‘ ) Sll- © 



Employing stress-strain (4.11) and integrating with respect to r 

E Uy ” foi+VMcxjr v!o 1 y 2 + a' 1 V :2 ]-^[3\D ( v- 1 - a' v v' 2 J + 

+|V ^°r] ^ cos © * fc©) (B.ll) 

Employing (4.12) and integrating with respect to 0 

Eu e 1 t" (^, + '&',) lo^v- * * *1 r 3 , 4 - lo/ ]sU .9 * 

- V>^- Vs/ lo^v + (a, +'o x 1 ')- k, r **] 

- + fCO (B.12) 



Prom (Boll) 



JL 

v 5© 



c> U& _ 

©V- 



- [(«!+ k.'V'?'’ + fc,v + a', v" 5 ]su.e + , 

'T’lw-e * h -p 

(B.13) 

[- + lot,- -2a',<-‘ , ]s^a + 

*» U,v* 2t,v- - 2aIf- 5 3 su.e *P' 

(B.14) 



Prom (B.12) 
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And also from (B.12) 

'f "7 s - l’ 0*1+ b, 1 ') “Sb.tr + <\ V 3 + Id, y '] «,u*. © 

+ P ^ - t, 1 +• (<*>+ t, 1 ) v- ' - b>,r - <N v * \ S'** ® 

+ f KQ)d9 _ _F (B.15) 

) y y. 

Equating (4.8) and (4.5) and using (B.13), (B.l4), and (B.15) 

+ [4t: A (!-«)«,] s-s/v 

or •$'(©) » fac®)4© <■ (vC'-f) a [4k,' »(>-«)&,] •5U.0 

(B.16) 

(vF-F l ) must necessarily be zero 

vF 1 - F = O => F=Cr 
Differentiating with respect to © 

5 '* 4 * [^b/ 4 (l-P^ft^CoS © (B.17) 

The solution to the differential equation (B.17) 

£ (.©) 3 A su. Q + (?>e©s © + 1.2 Is/ + ^ a ' ] ® sC Q 



Wavenumber > 2 
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From [29] 



. / Hi I -v\ . I -H+l\ 

\ - \<X*r + <o H \r ) to^e 



C v * - Iv%(^"l)a. w >r n V (w a -v> >l)V >H r V ' + vi(\a+i)(x'v,v ■*• 

+ ( W 1 +■ v\ - 2) b* v’* 3 Cos w 0 

6© * tv\(vA-0ct K v >u ’\ h(w *i) a! v" 1 *’ 1 * 

+ (h-i^(ia-x) b vi y " ^ cos h © 

6 V ©- v v\(vA+l^!ov,y - v*(vwl)ai r 

- V\(v\-l') lov* V 

Again employing (4.11) and integrating 

E u v * - [*>*.*""♦ c»-i) *••“ t u nK ;" “] • 

• Cos v* Q - vi [ w a.* r 1 '"’ + (u+i) bn v MA iAO. 1 , v ’*" 1 

- by» r~'] Cos u © +• 4 C©) 



Employing (4.12) and integrating 



(B.18) 



Eu© 3 ^v\o. u r H '+ (vi*4)l 6 *v**\v»a l ,v’ V * V (v»-4)lol V-^ 

, ^ r K-t 1 . \*4-l V -H -1 

• SUa u © 4- W V* + U D W V" -Ha Ct* V 



From (B.l8) JL 



£ ^ Uv 
V- d © 



From (B.19) £ 



+ *bLv" M+ '] scuu© -(Us)ae+Kv) 

(B. 19 ) 

C V" + Vv(w-l)b K v- V '-w a <Xu r U a - vv(vwO by, * 

. v' U ] 4 -V 5 V ' U " lj - V't U+ ‘ i ) bv>v"- rfo* V ^ 

- vn(u -* 0 C v“* 3^ S«l* ^ + 5- 

^ (bT 20 ) 

■ C[*C^K* + (wa-iX^+4) Id** 1 - vi(u*0 A* v” V '" 7 ’ 

- (v\-iYu-4)Io'k V-' W J 4-P •hw(tA-Vl) bvi ‘ 

• r'* — Vv £ U4 f) v* - Vv c^“0 ^o'w ^ J^S^aaK© *■ P ' 

(Bo 21 ) 
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Also from (B.19) 

0© _ 
"c v a 



4- f V\A* Y* + 

^ r . . h-^ 1 w [ -n-i 

.%U*w© 4. W> L VvCt ‘' Y 4-Vvb u V“ 4 -V\<X U V- 

4 v> W K V '*]$***© * ^ ^4C€>U& -£ (B.22) 



Equating (4.8) and (4.9) and using (B.20), (B.21), & (B.22) 



4 - F * + ^ - — ® O 

y 

^ = F-vY = cok^tai^t 

F * C 4- Or 

r A 5 i»* © 4 * B 6 £>T> © 



These general solutions must "be employed with the correct 
"boundary conditions to find the unique constants associated 
with each wavenumber. 
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APPENDIX C 

Buckling of a Circular Ring Subject to Stresses Uniformly 
Distributed Around the Edge 

The loading applied to the condition depicted in Figure 
C-l is analyzed for the first buckling mode. Assuming the 
deflection surface is a surface of revolution, Timoshenko 
[28] gives the required differential equation: 





i d <t> j _ Qr 2, 

r dv l dr Y d 


(0.1) 


where 


Q is the shearing force/length 




and D 


is the flexural rigidity 






Q = Nfc aw <}> = 


(0.2) 


Let 


N* 2 

K 3 ex! 


( 0 . 3 ) 




D 


(C. 3 ) 


becomes 








(C.K) 


If 


u.-= otv- (variable change) 




Then 


1 cL'b at & / 2 ^ 

U. + u 7£ + (u “ l ''° 


( 0 . 5 ) 



The general solution of which is 



+ ■ A,T,(u> AiY,M (0.6) 

Applying boundary conditions 

<^(v«o.)= o clamped inner edge (C.7) 

( sli 4. p —V O v-a b free outer edge ( C * 8 ) 

\4v + J 
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O 




Figure C-l 

Buckling of a thin circular plate around a rigid 
inclusion 



l6o 



Since 


^7 ■ A,o< T, (u) + Aj> Y ( *0 


(0.9) 


The boundary conditions (C.7)* (C.8) become 






A, It (*&•) + » 0 


(C.10) 




-»* ( A, + AiY l '(»<'»)) + 

+ -^(A,T,(^b) +AiY,(*'0): 0 


(C.ll) 


(C.10) is 


A .-a 


(0.12) 


And (C.ll) 


Y, ' 1 


(0.13) 


■JZ. T| = ^ T ’ ( 0 . 1*0 

(C.12), (C.13), (C.l4) are combined into a single equation. 
The first zero appears at about a value of .42 


To apply 


a 1 ^ » 
°- 3 * > b = 20 » ^ ' 3 




Where 


*-■41/0. then N c *".n&D 






£} = ^ 3 the flexural rieiditv 

12(1-*) 




Yields 


(5 cr s — = .oifcS EUV<x l 





h 
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